## figure 1 by subgroup

library(haven)
library(tidyverse)
library(weights)
library(labelled)
library(ggpubr)

rm(list=ls())
gss <- read_rds("clean-GSS-vars-w-weights.RDS")


#### PARTY


# select vars of interest and compute weighted mean for finrela and satfin 
gss1 <- gss %>% 
  dplyr::select(year, id, wtssps, satfin, finalter, finrela, rincblls, partyid, race,
                sex) %>%
  remove_labels() %>%
  # subset(is.na(mode) | mode==1)%>%
  mutate(satfin = as.numeric(dplyr::recode(satfin, 
                                    "1" = "3",
                                    "2" = "2",
                                    "3" = "1")))%>%
  
  mutate(pid3 = as.character(dplyr::recode(partyid, 
                                    "0"="Democrat",
                                  "1" = "Democrat",
                                    "2" = "Democrat",
                                    "3" = "Independent",
                                  "4"="Republican",
                                  "5"="Republican",
                                  "6"="Republican",
                                  "7"="other")))%>%
  subset(pid3 %in% c("Democrat", "Republican")) %>%
  group_by(year,partyid) %>%
  filter(is.na(finrela) == F & is.na(satfin) == F) %>% 
  dplyr::mutate(wtavg_finrela = weighted.mean(finrela, wtssps, na.rm=T), 
                lwr_finrela = wtavg_finrela-1.96*sd(finrela*wtssps, na.rm=T)/sqrt(n()), 
                upr_finrela = wtavg_finrela+1.96*sd(finrela*wtssps, na.rm=T)/sqrt(n()),
                wtavg_satfin = weighted.mean(satfin, wtssps, na.rm=T), 
                lwr_satfin = wtavg_satfin-1.96*sd(satfin*wtssps, na.rm=T)/sqrt(n()), 
                upr_satfin = wtavg_satfin+1.96*sd(satfin*wtssps, na.rm=T)/sqrt(n()))

gss_finrela <- gss1 %>% select(year, wtavg_finrela, lwr_finrela, upr_finrela,pid3) 
gss_finrela <- gss_finrela[!duplicated(gss_finrela), ]

gss_satfin <- gss1 %>% select(year, wtavg_satfin, lwr_satfin, upr_satfin,pid3) 
gss_satfin <- gss_satfin[!duplicated(gss_satfin), ]

finrela <- gss1 %>%
  select(finrela) %>%
  drop_na()

## Plot weighted average finrela 
a<- ggplot(gss_finrela, aes(x=year, y=wtavg_finrela, color=pid3)) +
  geom_point() +
  geom_line(linetype="dashed") +
  scale_y_continuous(limits = c(1,5),
                     labels = c("Far Below Average", "Below Average", "Average",
                                "Above Average", "Far Above Average"))+
  theme_bw()+
  scale_color_manual(values = c("grey", "grey3"))+
  labs(y = "Category", x = "", 
       title = "A", color = "Party",
       subtitle = "Compared with American families in general, would you say 
your family income is far below average, below average, 
average, above average, or far above average? N=70,428 (GSS)
Mean and 95% Confidence Intervals") +  
  geom_errorbar(aes(ymin = lwr_finrela, ymax = upr_finrela), linewidth=0.4) 
a


b<- ggplot(gss_satfin, aes(x=year, y=wtavg_satfin, color=pid3)) +
  geom_point() +
  geom_line(linetype="dashed") +
  theme_bw()+
  scale_y_continuous(limits = c(1,3), breaks=c(1,2,3),
                     labels = c("Not satisfied at all",
                                "More or less satisfied",
                                "Pretty well satisfied"))+
  scale_color_manual(values = c("grey", "grey3"))+
  labs(y = "Category", x = "", color="Party",
       title = "B", 
       subtitle = "Would you say that you are pretty well satisfied with your
present financial situation, more or less satisfied, or 
not satisfied at all? N=70,428 (GSS)
Mean and 95% Confidence Intervals") +  
  geom_errorbar(aes(ymin = lwr_satfin, ymax = upr_satfin), linewidth=0.4)
b

health_df <- gss %>%
  mutate(pid3 = as.character(dplyr::recode(partyid, 
                                           "0"="Democrat",
                                           "1" = "Democrat",
                                           "2" = "Democrat",
                                           "3" = "Independent",
                                           "4"="Republican",
                                           "5"="Republican",
                                           "6"="Republican",
                                           "7"="other")))%>%
  subset(pid3 %in% c("Democrat", "Republican")) %>%
  mutate(health = as.numeric(dplyr::recode(health, 
                                    "1" = "4",
                                    "2" = "3",
                                    "3" = "2",
                                    "4" = "1"))) %>%
  mutate(health_fac = factor(dplyr::recode(health, 
                                    "4" = "Excellent",
                                    "3" = "Good",
                                    "2" = "Fair",
                                    "1" = "Poor"),
                             levels = c("Excellent",
                                        "Good",
                                        "Fair",
                                        "Poor")))%>%
  mutate(happy = as.numeric(dplyr::recode(happy, 
                                   "1" = "3",
                                   "2" = "2",
                                   "3" = "1"))) %>%
  mutate(happy_fac = factor(dplyr::recode(happy, 
                                   "3" = "Very happy",
                                   "2" = "Pretty happy",
                                   "1" = "Not too happy"),
                            levels = c("Not too happy",
                                       "Pretty happy",
                                       "Very happy")))

health_df1 <- health_df %>%  
  group_by(year, pid3) %>%
  summarise(happy_w_av_yr = weighted.mean(happy, wtssps, na.rm=TRUE),
            sd_happy = sd(happy, na.rm=TRUE),
            se_happy = sd_happy / sqrt(n()),
            health_w_av_yr = weighted.mean(health, wtssps, na.rm=TRUE),
            sd_health = sd(health, na.rm=TRUE),
            se_health = sd_health / sqrt(n()),
            health_fac=health_fac) %>% unique()


c<- ggplot(health_df1, aes(x=year, y=health_w_av_yr, color=pid3)) + geom_point() +
  geom_line(linetype="dashed") +
  theme_bw() +
  scale_color_manual(values = c("grey", "grey3"))+
  labs(x="", y="Category",title="C", color="Party",
       subtitle="Would you say your own health, in general, 
is excellent, good, fair, or poor? N=58,448 (GSS)
Mean and 95% Confidence Intervals") +
  geom_errorbar(aes(ymin = (health_w_av_yr -1.96*se_health), 
                    ymax = (health_w_av_yr + 1.96*se_health)), width = 0.4) +
  scale_y_continuous(limits = c(1,4),
                     labels = c("Poor", "Fair", "Good", "Excellent"))
c


d<- ggplot(health_df1, aes(x=year, y=happy_w_av_yr, color=pid3)) + geom_point() +
  geom_line(linetype="dashed") +
  theme_bw() +
  scale_color_manual(values = c("grey", "grey3"))+
  labs(x="", y="Category",title="D", color="Party",
       subtitle="Taken all together, how would you say things are these days:
would you say that you are very happy, pretty happy,
or not too happy? N=70,869 (GSS)
Mean and 95% Confidence Intervals") +
  geom_errorbar(aes(ymin = (happy_w_av_yr -1.96*se_happy), 
                    ymax = (happy_w_av_yr + 1.96*se_happy)), width = 0.4) +
  scale_y_continuous(limits = c(1,3), breaks=c(1,2,3),
                     labels = c("Not too happy",
                                "Pretty happy",
                                "Very happy"))
d


png(file = "PNAS/figure-1-panel-by-party.png",
    height = 8.5,
    width = 14,
    units = "in",       # Use inches to match pdf dimensions
    res = 300)          # Set resolution (DPI) for high-quality output

print(ggarrange(a,b,c,d,
                nrow = 2, ncol = 2, common.legend = T, legend = "bottom"))  

dev.off()









##### BY SEX
gss2 <- gss1 %>%
  mutate(sex = as.character(dplyr::recode(sex, 
                                           `1` = "Male",
                                           `2` = "Female")))%>%
  group_by(year,sex) %>%
  filter(is.na(finrela) == F & is.na(satfin) == F) %>% 
  dplyr::mutate(wtavg_finrela = weighted.mean(finrela, wtssps, na.rm=T), 
                lwr_finrela = wtavg_finrela-1.96*sd(finrela*wtssps, na.rm=T)/sqrt(n()), 
                upr_finrela = wtavg_finrela+1.96*sd(finrela*wtssps, na.rm=T)/sqrt(n()),
                wtavg_satfin = weighted.mean(satfin, wtssps, na.rm=T), 
                lwr_satfin = wtavg_satfin-1.96*sd(satfin*wtssps, na.rm=T)/sqrt(n()), 
                upr_satfin = wtavg_satfin+1.96*sd(satfin*wtssps, na.rm=T)/sqrt(n()))

gss_finrela <- gss2 %>% select(year, wtavg_finrela, lwr_finrela, upr_finrela,sex) 
gss_finrela <- gss_finrela[!duplicated(gss_finrela), ] %>%
  drop_na()

gss_satfin <- gss2 %>% select(year, wtavg_satfin, lwr_satfin, upr_satfin,sex) 
gss_satfin <- gss_satfin[!duplicated(gss_satfin), ]%>%
  drop_na()

a<- ggplot(gss_finrela, aes(x=year, y=wtavg_finrela, color=sex)) +
  geom_point() +
  geom_line(linetype="dashed") +
  scale_y_continuous(limits = c(1,5),
                     labels = c("Far Below Average", "Below Average", "Average",
                                "Above Average", "Far Above Average"))+
  theme_bw()+
  scale_color_manual(values = c("grey", "grey3"))+
  labs(y = "Category", x = "", 
       title = "A", color = "Sex",
       subtitle = "Compared with American families in general, would you say 
your family income is far below average, below average, 
average, above average, or far above average? N=70,428 (GSS)
Mean and 95% Confidence Intervals") +  
  geom_errorbar(aes(ymin = lwr_finrela, ymax = upr_finrela), linewidth=0.4) 
a


b<- ggplot(gss_satfin, aes(x=year, y=wtavg_satfin, color=sex)) +
  geom_point() +
  geom_line(linetype="dashed") +
  theme_bw()+
  scale_y_continuous(limits = c(1,3), breaks=c(1,2,3),
                     labels = c("Not satisfied at all",
                                "More or less satisfied",
                                "Pretty well satisfied"))+
  scale_color_manual(values = c("grey", "grey3"))+
  labs(y = "Category", x = "", color="Sex",
       title = "B", 
       subtitle = "Would you say that you are pretty well satisfied with your
present financial situation, more or less satisfied, or 
not satisfied at all? N=70,428 (GSS)
Mean and 95% Confidence Intervals") +  
  geom_errorbar(aes(ymin = lwr_satfin, ymax = upr_satfin), linewidth=0.4)
b

health_df <- gss %>%
  mutate(health = as.numeric(dplyr::recode(health, 
                                           "1" = "4",
                                           "2" = "3",
                                           "3" = "2",
                                           "4" = "1"))) %>%
  mutate(health_fac = factor(dplyr::recode(health, 
                                           "4" = "Excellent",
                                           "3" = "Good",
                                           "2" = "Fair",
                                           "1" = "Poor"),
                             levels = c("Excellent",
                                        "Good",
                                        "Fair",
                                        "Poor")))%>%
  mutate(happy = as.numeric(dplyr::recode(happy, 
                                          "1" = "3",
                                          "2" = "2",
                                          "3" = "1"))) %>%
  mutate(happy_fac = factor(dplyr::recode(happy, 
                                          "3" = "Very happy",
                                          "2" = "Pretty happy",
                                          "1" = "Not too happy"),
                            levels = c("Not too happy",
                                       "Pretty happy",
                                       "Very happy")))%>%
mutate(sex = as.character(dplyr::recode(sex, 
                                        `1` = "Male",
                                        `2` = "Female")))
  
health_df1 <- health_df %>%  
  group_by(year, sex) %>%
  summarise(happy_w_av_yr = weighted.mean(happy, wtssps, na.rm=TRUE),
            sd_happy = sd(happy, na.rm=TRUE),
            se_happy = sd_happy / sqrt(n()),
            health_w_av_yr = weighted.mean(health, wtssps, na.rm=TRUE),
            sd_health = sd(health, na.rm=TRUE),
            se_health = sd_health / sqrt(n()),
            health_fac=health_fac) %>% unique() %>%
  drop_na()


c<- ggplot(health_df1, aes(x=year, y=health_w_av_yr, color=sex)) + geom_point() +
  geom_line(linetype="dashed") +
  theme_bw() +
  scale_color_manual(values = c("grey", "grey3"))+
  labs(x="", y="Category",title="C", color="Sex",
       subtitle="Would you say your own health, in general, 
is excellent, good, fair, or poor? N=58,448 (GSS)
Mean and 95% Confidence Intervals") +
  geom_errorbar(aes(ymin = (health_w_av_yr -1.96*se_health), 
                    ymax = (health_w_av_yr + 1.96*se_health)), width = 0.4) +
  scale_y_continuous(limits = c(1,4),
                     labels = c("Poor", "Fair", "Good", "Excellent"))
c


d<- ggplot(health_df1, aes(x=year, y=happy_w_av_yr, color=sex)) + geom_point() +
  geom_line(linetype="dashed") +
  theme_bw() +
  scale_color_manual(values = c("grey", "grey3"))+
  labs(x="", y="Category",title="D", color="Sex",
       subtitle="Taken all together, how would you say things are these days:
would you say that you are very happy, pretty happy,
or not too happy? N=70,869 (GSS)
Mean and 95% Confidence Intervals") +
  geom_errorbar(aes(ymin = (happy_w_av_yr -1.96*se_happy), 
                    ymax = (happy_w_av_yr + 1.96*se_happy)), width = 0.4) +
  scale_y_continuous(limits = c(1,3), breaks=c(1,2,3),
                     labels = c("Not too happy",
                                "Pretty happy",
                                "Very happy"))
d

png(file = "PNAS/figure-1-panel-by-sex.png",
    height = 8.5,
    width = 14,
    units = "in",       # Use inches to match pdf dimensions
    res = 300)          # Set resolution (DPI) for high-quality output

print(ggarrange(a,b,c,d,
                nrow = 2, ncol = 2, common.legend = T, legend = "bottom"))  

dev.off()



######### WHITE/BLACK

gss3 <- gss %>% 
  dplyr::select(year, id, wtssps, satfin, finalter, finrela, rincblls, partyid, race,
                sex) %>%
  remove_labels() %>%
  # subset(is.na(mode) | mode==1)%>%
  mutate(satfin = as.numeric(dplyr::recode(satfin, 
                                           "1" = "3",
                                           "2" = "2",
                                           "3" = "1")))%>%
  
  mutate(race = as.character(dplyr::recode(race, 
                                           "1" = "White",
                                           "2" = "Black",
                                           "3" = "Other")))%>%
  subset(race %in% c("White", "Black")) 

df <- gss3%>%
  group_by(year,race) %>%
  filter(is.na(finrela) == F & is.na(satfin) == F) %>% 
  dplyr::mutate(wtavg_finrela = weighted.mean(finrela, wtssps, na.rm=T), 
                lwr_finrela = wtavg_finrela-1.96*sd(finrela*wtssps, na.rm=T)/sqrt(n()), 
                upr_finrela = wtavg_finrela+1.96*sd(finrela*wtssps, na.rm=T)/sqrt(n()),
                wtavg_satfin = weighted.mean(satfin, wtssps, na.rm=T), 
                lwr_satfin = wtavg_satfin-1.96*sd(satfin*wtssps, na.rm=T)/sqrt(n()), 
                upr_satfin = wtavg_satfin+1.96*sd(satfin*wtssps, na.rm=T)/sqrt(n()))

gss_finrela <- df %>% select(year, wtavg_finrela, lwr_finrela, upr_finrela,race) 
gss_finrela <- gss_finrela[!duplicated(gss_finrela), ]

gss_satfin <- df %>% select(year, wtavg_satfin, lwr_satfin, upr_satfin,race) 
gss_satfin <- gss_satfin[!duplicated(gss_satfin), ]


a<- ggplot(gss_finrela, aes(x=year, y=wtavg_finrela, color=race)) +
  geom_point() +
  geom_line(linetype="dashed") +
  scale_y_continuous(limits = c(1,5),
                     labels = c("Far Below Average", "Below Average", "Average",
                                "Above Average", "Far Above Average"))+
  theme_bw()+
  scale_color_manual(values = c("grey", "grey3"))+
  labs(y = "Category", x = "", 
       title = "A", color = "Race",
       subtitle = "Compared with American families in general, would you say 
your family income is far below average, below average, 
average, above average, or far above average? N=70,428 (GSS)
Mean and 95% Confidence Intervals") +  
  geom_errorbar(aes(ymin = lwr_finrela, ymax = upr_finrela), linewidth=0.4) 
a


b<- ggplot(gss_satfin, aes(x=year, y=wtavg_satfin, color=race)) +
  geom_point() +
  geom_line(linetype="dashed") +
  theme_bw()+
  scale_y_continuous(limits = c(1,3), breaks=c(1,2,3),
                     labels = c("Not satisfied at all",
                                "More or less satisfied",
                                "Pretty well satisfied"))+
  scale_color_manual(values = c("grey", "grey3"))+
  labs(y = "Category", x = "", color="Race",
       title = "B", 
       subtitle = "Would you say that you are pretty well satisfied with your
present financial situation, more or less satisfied, or 
not satisfied at all? N=70,428 (GSS)
Mean and 95% Confidence Intervals") +  
  geom_errorbar(aes(ymin = lwr_satfin, ymax = upr_satfin), linewidth=0.4)
b



health_df <- gss %>% mutate(race = as.character(dplyr::recode(race, 
                                                               "1" = "White",
                                                               "2" = "Black",
                                                               "3" = "Other")))%>%
  subset(race %in% c("White", "Black"))%>%
  mutate(health = as.numeric(dplyr::recode(health, 
                                           "1" = "4",
                                           "2" = "3",
                                           "3" = "2",
                                           "4" = "1"))) %>%
  mutate(health_fac = factor(dplyr::recode(health, 
                                           "4" = "Excellent",
                                           "3" = "Good",
                                           "2" = "Fair",
                                           "1" = "Poor"),
                             levels = c("Excellent",
                                        "Good",
                                        "Fair",
                                        "Poor")))%>%
  mutate(happy = as.numeric(dplyr::recode(happy, 
                                          "1" = "3",
                                          "2" = "2",
                                          "3" = "1"))) %>%
  mutate(happy_fac = factor(dplyr::recode(happy, 
                                          "3" = "Very happy",
                                          "2" = "Pretty happy",
                                          "1" = "Not too happy"),
                            levels = c("Not too happy",
                                       "Pretty happy",
                                       "Very happy")))

health_df1 <- health_df %>%  
  group_by(year, race) %>%
  summarise(happy_w_av_yr = weighted.mean(happy, wtssps, na.rm=TRUE),
            sd_happy = sd(happy, na.rm=TRUE),
            se_happy = sd_happy / sqrt(n()),
            health_w_av_yr = weighted.mean(health, wtssps, na.rm=TRUE),
            sd_health = sd(health, na.rm=TRUE),
            se_health = sd_health / sqrt(n()),
            health_fac=health_fac) %>% unique()


c<- ggplot(health_df1, aes(x=year, y=health_w_av_yr, color=race)) + geom_point() +
  geom_line(linetype="dashed") +
  theme_bw() +
  scale_color_manual(values = c("grey", "grey3"))+
  labs(x="", y="Category",title="C", color="Race",
       subtitle="Would you say your own health, in general, 
is excellent, good, fair, or poor? N=58,448 (GSS)
Mean and 95% Confidence Intervals") +
  geom_errorbar(aes(ymin = (health_w_av_yr -1.96*se_health), 
                    ymax = (health_w_av_yr + 1.96*se_health)), width = 0.4) +
  scale_y_continuous(limits = c(1,4),
                     labels = c("Poor", "Fair", "Good", "Excellent"))
c


d<- ggplot(health_df1, aes(x=year, y=happy_w_av_yr, color=race)) + geom_point() +
  geom_line(linetype="dashed") +
  theme_bw() +
  scale_color_manual(values = c("grey", "grey3"))+
  labs(x="", y="Category",title="D", color="Race",
       subtitle="Taken all together, how would you say things are these days:
would you say that you are very happy, pretty happy,
or not too happy? N=70,869 (GSS)
Mean and 95% Confidence Intervals") +
  geom_errorbar(aes(ymin = (happy_w_av_yr -1.96*se_happy), 
                    ymax = (happy_w_av_yr + 1.96*se_happy)), width = 0.4) +
  scale_y_continuous(limits = c(1,3), breaks=c(1,2,3),
                     labels = c("Not too happy",
                                "Pretty happy",
                                "Very happy"))
d


png(file = "PNAS/figure-1-panel-by-race.png",
    height = 8.5,
    width = 14,
    units = "in",       # Use inches to match pdf dimensions
    res = 300)          # Set resolution (DPI) for high-quality output

print(ggarrange(a,b,c,d,
                nrow = 2, ncol = 2, common.legend = T, legend = "bottom"))  

dev.off()




